STABLE NUMERICAL METHODS FOR TWO CLASSES OF SDEs 
WITH MULTIPLICATIVE NOISE: BILINEAR AND SCALAR 



H. A. MARDONES *+ AND C. M. MORA ** 

Abstract. We develop new numerical schemes for bilinear systems of stochastic differential 
equations (SDEs). To this end, we present a new methodology for solving stiff multidimensional 
SDEs with multiplicative noise, as well as, we introduce a new stable numerical method for non-linear 
scalar SDEs. The rate of weak convergence of the new schemes is linear. Moreover, they preserve the 
possible exponential stability of the unknown solutions for any step-size. Four numerical experiments 
illustrate the good performance of the proposed algorithms. 
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1. Introduction. This paper introduces a new approach to solve stiff stochastic 
differential equations (SDEs) with multiplicative noise, namely, SDEs of the form 

pt m pt 

X t =X + b(X s )ds + y2 a k {X s )dW* (1.1) 
Jo k=1 Jo 

whose numerical solutions by the Eulcr-Maruyama scheme exhibit incorrect behav- 
iors. Here, W 1 , . . . , W m arc independent real valued Wiener processes on a filtered 
complete probability space ($t) t >o i^ji X t is an adapted Revalued stochas- 

tic process, b, a k : M. d — > M. d have continuous first-order partial derivatives, and 
Xq G L 2 (fl,P). For simplicity, we only consider autonomous SDEs. Our final goal is 
to have a set of schemes that preserve dynamical properties of (|1.1[) , like asymptotic 
stability, for any step size A of the time discretization. 

Section [2] is devoted to scalar SDEs. More precisely, we develop a new scheme 
for (fTTTjl with d= 1. In case 6(0) =0 and a k (0) = 0, the new numerical method 
preserves, for any step size A > 0, both the sign of Xq and the possible exponential 
stability of X t , that is, the property limsup^^ -|-ln||A" t || < —A < . This paves 
way for the main objective of this article: to provide stable schemes for computing 
the mean value of / (Xt), where / : M. d — > K is smooth and, by abuse of notation, 

pt rn nt 

X t =X Q + BX s ds + J2 ° k X s dW k , (1.2) 
Jo k=1 Jo 

with B, a 1 , . . . , er m real matrices of dimension d x d. In this research direction, Sub- 
section [XT] presents a general technique for constructing stable methods for (|1.1[) with 
d > 1. Then, in Subsections 11.21 and 13.21 we use the new ideas to derive promising 
weak schemes for (|1.2[) . The bilinear SDE (|1.2|) arises, for example, from the spa- 
tial discretization of some stochastic partial differential equations (see, e.g., [T| IT2"]). 
and describes important dynamical features of non-linear SDEs via the linearization 
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around their equilibrium points (see, e.g., 0133]). Furthermore, the new numerical 
methods for (|1.2[) guide us in the application of our general methodology to systems 
of non-linear SDEs. 

1.1. Previous works. In many cases, the semi-implicit and explicit Euler meth- 
ods preserve dynamical properties of the underlying SDEs provided that the step size 
of the discretization is small enough. Recall that (see, e.g., [T71 119] ) 

limsup - In ||A t j < — A a.s. (1.3) 

t— >oo t 



whenever 



||6(x)|| <K n \\x\\ V||z|| <nandV?ieN, 

||cr fe (x)\\ < K\\x\\ VxGlR d , and (1.4) 

_ A;= mp f<*H*)) + iTZ, . !%.,(»■<*<*»") <ft (1 . 5) 

xei d ,i/o \ ||.r|| |.t| J 

Here and subsequently, K stands for a generic positive constant and K n > 0. If 
Conditions (|1.4[) and (|1.5|l hold, together with b (0) = 0, then Higham, Mao and Yuan 
[17] proved that the Euler-Maruyama scheme 

m 

E n+1 =E n + b (E n ) A + (En) (V(t+1)A ~ ^nA 

fe=l 

is almost sure exponentially stable for sufficiently small step sizes A > in case 

||6(ic)|| < ir||a;|| VxeM d . 
At the cost of solving systems of algebraic equations, the backward Euler method 

m 

E n+1 =E n + b (E n+1 ) A + Y,*" (En) (W ( fc „ +1)A - WZ A ) (1.6) 

k=l 

improves the numerical stability of E n (see, e.g., [TC] [TTJ [35] ) . For instance, Higham, 
Mao and Yuan [T7] established that there exists Ao > such that for any A € ]0, Ao], 
E n is almost sure exponentially stable provided that (|1.4[l is valid, b (0) = 0, and 

(,-,,K,)- tw > +8UP fSkKwi! _ iz^ (*))') <0 . (L7) 

J, \\x-yf 2||i-|| 2 H y 

The Euler schemes i?„ and £"„ exhibit a poor numerical performance in situations 
where, for example, some partial derivatives of the diffusion coefficients a k are not 
small. A simple model problem for such SDEs is 

dX t = XX t dW t , (1.8) 

with A > (see, e.g., [121 (23]); the trajectories of E n and E n applied to (|1.8[) blow 
up unless A is very small. Using (|1.8[) as a motivational problem, Milstcin, Platen 
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and Schurz |22] introduced the general formulation of the balanced implicit methods, 
a class of fully implicit schemes for whose implementation depends of the choice 
of certain weights (see, e.g., [H [T5J [H]). To the best of our knowledge, the reported 
balanced schemes have good asymptotic stability properties and low speed of weak 
convergence. 

Other implicit integrators for together with their predictor-corrector ver- 

sions, arise from the Ito- Taylor expansions of Xt (see, e.g., [13j [23j [27l [2SJ |30l 136]). 
In particular, Kloeden and Platen [T5] proposed the following class of weak implicit 
schemes: 



E n+1 = E n + (a c [E n+1 j + (1 - a) c [E n j J A (1.9) 

m 

+ (V (En+l) + (1 - V) ° k (En)) VZ£ 



k=l 



where q,i] £ [0,1], : k = l,...,mandn > 0} is a collection of independent 
random variables taking values ±1 with probability 1/2, and 

m d da k 

C(X) =b(x)- V Y,J2' TJ:k ^faT^- 

We have E n rj X ni \- Applying (|1.9[) with a = r\ = 1 to (| 1 . 8[) yields a fully implicit 
method which is almost sure asymptotically stable, but converges to as n — > oo 
too much faster than X t . Taking a = r\ = 1/2 gives a trapezoidal scheme with good 
asymptotically stability properties, nevertheless it fails to preserve the sign of Xq in 
the numerical solution of (|1.8p . 

Finally, numerical methods adapted to specific types of SDEs with multiplicative 
noise have been developed, for instance, in [3j [23l EH |25] . The numerical integration 
of mean-square stable SDEs has been treated, for example, in [TJ [TU1 QTJ HS1 EH] • To the 
best of our knowledge, this paper is the first work to present numerical methods for 
SDEs that preserve, for any step size A > 0, the asymptotic stability of the solutions 
of relevant classes of SDEs (see, e.g., P^lfTT] ). 

1.2. New numerical method for bilinear SDEs. It follows from 

x < = Wi ml 

that we can divide the computation of the solution of (|1.2|) into: 

(i) The approximation of X t := X t j \ X t \ by an adapted stochastic process tak- 
ing values in the unit sphere. 

(ii) The numerical simulation of \X t \ by means of a scheme that preserves the 
dynamical properties of ||X t ||. 

Since f| 1 . 2 j) is bilinear, using Ito's formula yields 

X t =X + B (! s ) X s ds + Y, (° k - (X» ° k X s )) X s dW* (1.10) 



a k x f 



(see Subsection 13. II for details), where, by abuse of notation, 

B(X.) =B-(X S ,BX S )+Y^ U(x s ,a k X s ) -(x s ,a k X s )a k -- 

k=l 
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Although (|1.10j) is a locally Lipschitz SDE, its solution has norm 1 for all t > 0, 
and hence we can improve the performance of the numerical schemes applied to (|1.10[) 
by projecting on the unit sphere at each discretization step. This projection procedure 
has been used with success in the numerical solution of the non-linear Schodingcr 
equations (see, e.g., [24"I[26] ). 

Moreover, we utilize Ito's formula to obtain 

l|XtlM|Xoll+ io { ra *h \\ Xs f ) ds 

+ V f {Xs > akXs) dw k (in) 

h'o \\ x >\\ 

(see Subsection 13.11 for details). In Section [5] we develop a stable scheme for scalar 
stiff SDEs that shows very good performance in numerical experiments. A close look 
at this numerical method leads us to divide and multiply each integrand of (jl.lip by 
||X s j|, and so (jl.lll) becomes 

\\x t \\ = \\ x o\\+ 1 \ i x s' Bx °)+lit (lk feJ ^f ~ (*..***.> 2 )) ra^ 

m „ f 

+ V / (X s ,a k X s )\\X s \\dW k . (1.12) 
, Jo 



We propose to approximate numerically X t by solving the system formed by (| 1 . 10[) 
and (|1.12p . Next, we present a simple way to do this. For simplicity, from now on we 
consider the equidistant time discretization T n = nA, where A > and n = 0, 1, . . . 
Suppose that X n ps Xt„- Applying the weak Euler approximation to (|1.10j) we get 
Xr n+1 ~ Z n+ i with 



Z n+ i =X n + B (X n ) X n A + J2(' jk - ( X n,(r k X n )) X n VA 

fe=l 

where £o,£o, • • • i £o™j £i> • • • are independent random variables taking values ±1 with 
probability 1/2. Projecting Z n +i on the unit sphere gives 

X T n+1 ~ X n+ i := Z n+ i/ ||Z„ +1 || . (1-13) 

From (|1.12|) we deduce that ||At||, with t £ [T„,T n+ i], is well approximated by 
the solution of the linear scalar SDE 

T k = Vn + J T ( (X n ,BX n ) + -^(j\a k X4 3 -(X n ,a k X n ) 2 Uri a da (1.14) 

m „t 

+ X) / ( X n,a k X n ) Vs dW k , 



k=l 
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Fig. 1. Computation o/ E arctan ^1 + (X*) 2 ^ , where t £ [0,10] and Xt solves 11.161 1 with 
b = —2, cr = 4 and e = 4. The true values are plotted with a solid line. The circles (resp. stars) 
represent the approximations of Earctan ^1 + (X*) 2 ^ , with t = 0, 1, .. . , 10, given by f\ n X n (resp. 
the backward Euler method E n ). 

where fj n ps ||Xr n II ■ Replacing Wj. n+1 — W£ n by >/A£* in the explicit solution of (|1.14|) 
we get 

(/ j m m \ 

(X n ,BX n ) + -J2\\<r k X n \\ ~ ^(A > „,a fc A > „) 2 A 

rn \ 

+ £<X n)t x fe X n )v / A^ , (1.15) 
fc=i / 

and so ^„+i ps ||Xt„ +1 ||- Iterating (|1.13p and (|1 . 15[) gives the recursive scheme fj n X n . 

In Section[31 we prove that fj n X n approximates Xt k with weak rate of convergence 
equal to 1, and we establish that for any A > 0, 

limsup — — ln||?7„X n || < — A P — a.s. 

whenever (|1.2p satisfies the condition (|1.5p . Moreover, Section [3] develops versions of 
fj„X n adapted to systems where B is very ill-conditioned. 

We will illustrate the numerical behavior of fj n X n by means of the test equation 

X t = X + f (l f) X s ds + f ft °) X s dWl + /7° -f] X s dWl (1.16) 



x bj s J \0 a) s s J Q \e / 

which has been proposed in [HUH]- We take Xo = (1, 2) T , b = —2, a = 4 and e = 4, 
and so A = —2. Hence we are dealing with an exponentially stable SDE. To avoid 
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variance problems, we calcule Earctan ^1 + (A^ 1 ) 2 ^ , whose reference values (solid 
lines) have been obtained by sampling 10 s times the explicit solution of (|1.16[) . Here 

T / 2 

Xt = [Xt, X?) e R 2 . Figure Q] compares the computation of E arctan I 1 + (X^) 

by using fj n X„ (represented by circles) with that produced by the backward Euler 
method E„ defined by (|1.6p with W k n+1 ^ A — W k A replaced by V~At; k (represented by 
stars); all the sample sizes are equal to 10 6 . Figure [1] shows the very good qualitative 
behavior of fj n X n in the numerical solution of (|1.16p . We can observe that the first 
coordinate of fj n X n decays to with the same speed that the true solution, even 
for large step-sizes. In contrast, the trajectories of E n blows up when A = 0.1 and 
A = 0.02. Moreover, fj n X n achieves an excellent accuracy in cases A = 1 and A = 0.1. 

2. Stable schemes for scalar SDEs. In this section, we restrict our attention 
to stiff scalar SDEs, that is, we focus on 

pt m pt 

X t =X + / b(X s )ds + Y] / a k (X s )dW k , 
Jo fc=1 Jo 

where b, a k : R -s- R are continuously diffcrcntiable functions. 

2.1. Derivation of the numerical method. We begin by assuming 6(0) = 
cr 1 (0) = • ■ ■ = (J rn (0) = 0. Then 

x t = x + f b -^±x sds + jr f a -^x s dw k , 

Jo A s k=l Ja s 

where by abuse of notation, we write b (0) /0 and a k (0) /0 instead of the derivatives 
at of b and a k respectively. 

Suppose that X n is an 3t„ -measurable random variable such that X n sa Xt„, 
here and subsequently, Tj = jA for all j = 0, 1, . . . Then, for all t £ [T n , T n +i], 



Xt ~ X n + 



•J T n s J T n s 

Since cc i — ^ 6 (x) / x and x v-¥ a k (x) jx are continuous functions, (|2.1[) leads to 

ft b (X n ) ™ f* a k (X n ) 

Jt„ -X-n k=l Trl 

Hence X t is approximated by the solution of the linear scalar SDE 

r* b (X n ) ™ [* a k (X n ) , , 

Y t =X n + -\Hl Ys ds + y2 £ ' Y s dW k , (2.2) 

JT„ -X-n k=1 Jt„ X-n 

and so Xt„ +1 ~ Yt„ +1 ■ Solving explicitly (|2.2[) gives 



X„exp 
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As we are interested in weak approximations of X t , we replace Wj, — W£ 



in {33} by VAW*, where W \ Wq , . . . , W m , Wj 1 , . • • are independent identically dis- 
tributed (i.i.d.) random variables with symmetric law and variance 1. This gives the 
numerical scheme 

(2.4) 

From (|2.4I) it follows that X n preserves the sign of the initial data. Furthermore, 
we next establish that X n is almost sure exponential stable for any A > provided 
that er 1 , . . . , a m are globally Lipschitz and 

-A:= sup \b{x)/x- jr(a k (x)/x) 2 /2 ) < 0. 
x£S.,x?0 \ fe=1 y 

Hence X n preserves two important dynamical properties of the solution of (II. ip . 

Theorem 2.1. Consider \X n j n>Q given by the recursive formula ^2.4\ )- Suppose 

that \1.J$ and i fi.J)) hold, together with b (0) = and E (^o) 2 < °o- TTien 

lim sup — !— In |Xj < -A P - a.s. (2.5) 



A 



n 

Proof. Deferred to Section |2"XT1 □ 

The high performance achieved by X n in our numerical experiments, together 
with its good theoretical properties, motivates us to adapt X n to the framework 
where b (0) , a 1 (0) , . . . , a m (0) are not necessarily equal to 0. To this end, we rewrite 
(fl~T|) in the form 

X t = X + ^(^^X s + H 0)) d s 
m « 'a* (JT.) - a* (0) v , k 



+ E/ o { " \ W *. + »*(0)jdW?, 
and so for all i € [T„,T n+ i], 

X t « X„ + (n (X n ) X s + b (0)) ds + J2j T ( Xk X s + a k (0)) dW : 



k=l 

where X n is an -measurable random variable approximating Xt„, 

u^-f*-^, if^O Afe (x) (x ) - * k (0)) /x, if^O 

6'(0), ifx = andA K)'(0), ifx = ' 

This leads to approximate X t by the solution of 

pt m 

Y t =X n + ((j, (X n ) Y s + b (0)) ds + Y, (A fc (A > „) F a + <7 fe (0)) oW s fc . (2.6) 

« TU. J 1 J Tn 
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The explicit solution of (12.61) is 



(/ m \ ft rn „t 

X n +[b(0)-J2X k (X n )a k (0)) ^ds + ^Ho) K 1 
\ fc=i / Jt ™ k=i Jt " 



where 



$ t = exp ( (X n ) - 1 fj A fc (X„) 2 ^) (t - T n ) + fj A fc (X n ) (W k - W*J j . 
Since "Jj 1 « for all s <E [T„,T ra+ i], -X"t„ +1 is approximated by 

$ T „ +1 (x„ + (b (0) - fj A fc (X n ) a k (0)^ A + «r fc (0) (w* +1 - W k n ) J . 
This yields the weak scheme 
X n+1 = <S„ +1 ( X n + 1 6(0) - jr X k (X n ) a k (0) ) A + a k (0) VAW k ) , (2.7) 



fc=i / fc=i 



where Wq , Wq , . . . , Wo™, W"*, . . . are i.i.d. random variables with symmetric law and 
variance 1, and 

= CX P ^ L (Xn) ~\fl X " A + £ A " W ^) ' 



The following theorem establishes that X, given by (|2.7[) , converges weakly to X 
with order O (A) under a basic set of assumptions. In case b (0) = c 1 (0) = . . . = 
<7 m (0) = 0, (|2.7[) becomes (|2.4p . and so the rate of weak convergence of X, defined 
by (|2.4p . is equal to 1. 

Notation 2.1. We will use the same symbol K (•) (resp. K) for different positive 
increasing functions ( resp. positive real numbers ) having the common property to be 
independent of A. Similarly, q denotes generic constants greater than or equal to 
2. We write (R d ,R) for the set of all t-times continuously differ entiable functions 
f : R d — > R such that f and all its partial derivatives of orders 1,2, ... ,£ have at most 
polynomial growth. 

Theorem 2.2. Let b, a 1 , . . . , o~ m be Lipschitz continuous functions belonging to 

C^(R,R) such that \b(x)\ + jcr 1 (x)\ H h |er m (x)\ < K(l+ \x\) for all x G R. Fix 

T > and f G C£(K,R). Consider the scheme X described by with A = T/N , 

where N € N. Assume that Ecxp (rW k ^j < oo for all r > 0, Xq has finite moments 

of any order, and that for every g € Cp (R, R) , 

\Eg(X )-Eg(X )\ < K (1 + E |Jf | 9 ) T/N G N. 

T/ien /or a« iV e N, 

|E/ (X T ) - E/ (Xjv) | < K (T) (1 + E \X \ q ) T/N. (2.8) 
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Proof. Deferred to Section 12X21 □ 

Remark 2.1. IfW% are standard Normal random variables, then 
Eexp (rW^j = exp(r 2 /2) . 

In case W% are bounded random variables we also have Eexp ^rffjf j < oo. 

2.2. Numerical experiment. In this subsection, we illustrate the behavior of 
X n by means of the locally Lipschitz SDE 

X t = X + j [aX s - b (X,) 3 ) ds + J aX s dW}, (2.9) 

where b,o~ are real positive numbers, a € R, and Xq = 1. This scalar cubic SDE is 
known as the stochastic Ginzburg-Landau equation, and constitutes a classical test 
equation in the theory of stochastic bifurcation (see, e.g., [U [7]). Let £o;£o> ■ • ■ i£o\ 
£l, . . . be independent random variables taking the values ±1 each with probability 
1/2. Then, we will solve numerically (|2.9[) using three schemes: X n defined by (]2.4p 
with W% replaced by the backward Eulcr method E n given by (|1.6I) with vA^J 
in place of W[° n+1)A - W^ A , and 

K+Mi = K CX P ((« - ^ 2 /2)A + oyfEg) 

ZUx = (l - &A (^ +1/2 ) 2 /2) / (l + bA (^ +1/2 ) 2 /2) 



Fig. 2. Computation o/Eln ^1 + (Xt) 2 ^ , where t 6 [0, 5] and Xt satisfies 112.91 1 JuiWi a = 6 = 1 
and <t = 2. The "true" values are plotted with a solid line. The circles, stars and diamonds stand 
for the schemes X n , E n and Z s respectively. The step sizes 1, 0.5 and 0.25 are represented by 
dashdot, dashed and dotted lines respectively. 
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It is worth pointing out that E n entails the solution of a nonlinear equation at each 
step, and that is a weak version of the splitting-step algorithm for (|2.9[) introduced 
by Subsection 4.2 of f2"5] . 

Figures [21 [3] and Table POl show features of the computation of Eln ^1 + (X t ) 2 ^j 

obtained from the sample means of 10 s observations of X n , E n and Z^. The solid 
line identifies the "true" values gotten by sampling 10 8 times E n with A = 2~ n w 
0, 000488. The lengths of all the 99% confidence intervals arc at least of order 10 -3 , 
they have been estimated following [T5] , 

First, we take a = b = 1 and a = 2, which is the motivating example of [17] . 
Since b (x) /x - (a 1 (x) /xf /2 < -1, 

limsup - In \Xt\ < —1 a.s. 

t— >oo t 

From jTT] we have that the Euler-Maruyama scheme applied to (|2.9[) blows up, with 
positive probability, at a geometric rate. In numerical experiments, we saw that E n 
and fail to preserve the sign of Xq for A = 1 and A = 0.5. To this end, we 
computed 10E (tt/2 - arctan (l0 3 X t + 10 2 )). 

Figure [2] presents the numerical solution of (|2.9j) . with a = b = 1 and a = 2, using 
the step sizes A = 0.25,0.5, 1. Figure [5] suggests us that X n replicate very well the 
long time behavior of X t , even for A = 1. Moreover, we can see that the accuracy of 
X n is very good, even for large step sizes; X n achieves significantly lower errors than 
Z s , which is a method adapted to the characteristics of (|2.9p . 



Fig. 3. Computation of Eln ( 1 + (Xt) ) , where t e [0,10] and X t solves 1 12. 91 1 with a = 6, 
9 and a = 3. The "true" values are plotted with a solid line. The circles, stars and diamonds 



represent the schemes X„, E n and Z 3 respectively, 
dashed and dotted lines respectively. 



The step sizes 0.25 and 0.0625 are denoted by 



A = 0.25 




\ 

A = 0.25 



A = 0.0625 



0- 
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A 


|E/ (*r) -E/ (^t/a)I 


Ef (X T ) - Ef ( Zs.A 

./ \ i ) J y T/A / 


|E/ (X t ) -Ef(X T/A )\ 


1 


0.51106 


on nn'7'7 
29. OU/ 1 


0.0/461 / 


1/2 


0.51504 


26.1167 


0.05854 


1/4 


0.45996 


21.5922 


0.026918 


1/8 


0.0089311 


4.1708 


0.0078809 


2 -4 


0.048116 


0.068864 


0.003364 


2- 5 


0.015164 


0.041187 


0.0029212 


2 -6 


0.0056366 


0.022123 


0.0017523 


2" 7 


0.0023573 


0.011367 


0.00095535 



Table 2.1 

Estimation of errors involved in the computation offLf (Xt) forT = 10 and f (x) = In (l + x 2 ) . 
Here, Xt verifies 12. 91) with a = 6, 6 = 9, cr = 3 and Xq = 1. 



Second, we choose a = 6, 6 = 9 and a = 3. Then 

-A= sup (a- 6x 2 - ct 2 / 2 ) > 0, 

and so Condition (|1.5[) does not hold. In this case, (|2.9p has three invariant forward 
Markov measures (see, e.g., [J). Calculating 10E (jr/2 — arctan (l0 3 A" t + 10 2 )) we 
observe that E n (resp. can take negative values when A > 1/8 (resp. A > 1/16). 
Figure [3] displays the numerical approximation of Em ^1 + (X t ) 2 ^j by means of 

X n and E n with step sizes A = 0.25 and A = 0.0625, as well as by using with 
A = 0.0625. Moreover, Table HOI shows errors made in the weak numerical integration 

of (|2.9j) . where the reference value of E In ( 1 + (Xiq) 2 ) was obtained by sampling 10 s 



times E n with A = 2 11 . In this test problem, the new scheme X n again provides 
very good approximations of E/ (X t ), even for large values of A. 

2.3. Proofs. 

2.3.1. Proof of Theorem Hm From (3HJI it follows that 

H*H,i-Mi*,i + i:^-i£(^)V^- 




with S n = E?=o£fc=i ° k i X j) l X j vAW*. By (ESJ, we thus get 



— 5—lii|X n +i| < —J— lnlXol - AA + — |— 5„. (2.10) 
7i+l 1 1 n + 1 1 1 n + 1 

Using (TO} gives E f£™ =1 (-Xj) /^j v^IU/) 2 < ifA, and so 

\ 2 



V^_ E (V^B1VA^) <oo. 



Since 



fc=l j 
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applying a generalized law of large numbers we deduce that S n / (n + 1) — > a.s. (see, 
e.g., p. 243 of [H]). Then, letting n -> oo in ((2~TU|) we obtain (j23j) . □ 

2.3.2. Proof of Theorem 12.21 To shorten notation, for any n > we set 

/„ := I \i (X n 



1 m \ m 

fe=i / fe=i 



and 3 „ := (b (0) - £™ =1 A fc (X n ) a k (0)) A + £™ x a fc (0) VA#„ fc . Then 
X n+1 = cxp (/„) (l n + g n ) = X n + (e fn - 1 - f n ) X n + (e fn - l) g n + f n X n + g n , 
and hence 

n n n 

X n+1 =X + Y, (e /fc ~l-fk)X k + Y, ( efk - 1) 0* + E (A** + fffc) ■ 



fc=0 



fc=0 



Let <? > 2. By 



fc-i 

exp(x) - ^V/(j!) 

3=0 



< \x\ exp(|x|) , 



using Holder's inequality yields 



X n +i < 



K^ + Kin + lf- 1 [J2\M 2q ^ fkl \X k \ q + J2\fk\ q \9k\ q e^ 



\k=0 



-K (n + l)* -1 j^A 9 

k=Q 
n 

-K (n + l)* _1 $^A« 



1 m 



fe=0 



fe=0 



3=1 



XA? /2 £ (A J (X fe ) X, + a* (0)) W^' 



fe=0 



For any t > 0, 

Eexp ( 
and so for all I € N, 



< Eexp (tWi) +Eexp(-tW J k ) < oo, 



< i\ E cxp 



(N) 



< 00. 



(2.11) 



(2.12) 



(2.13) 



Since [i and A fc are bounded functions, we use (|2.12j) . (|2.13p and the Burkholder- 
Davis-Gundy inequality to obtain from (|2.11j) that 

n 

E \X n+1 \ q < KE | X \ q + K (T) + K (T) A ^ E \X k \ q , 

k=0 



STABLE SCHEMES FOR SDEs 



13 



with n = 0, . . . , N — 1. Applying a discrete Gronwall lemma (see, e.g., [I]) we deduce 
that for all n = 0, . . . , N , 

E\X n \ q <K(T)(l+E\X \ q ) . (2.14) 
Consider again q > 2. Using 

\X n+1 -X n \ < \e fn - l| \X n \ + e fn \g n \ < |/„| cxp(|/„|) \X n \ + exp(/„) \g n \ , 
together with (j2"7H2|) and ([2~T3jl . we get 

E(\X n+1 -X n \ q /$ Tn ) <K{T)A"/ 2 {l + \X n \ q ). (2.15) 

Here, we assume without loss of generality that ■ ■ -W™ are 3r„+i -measurable 
and independent of $t„ ■ 

Since \X n+l - (1 + /„ + /»/2 + (X„ + fln ) | < |/ n | 4 exp (|/ n |) |X n + <?„| , 

m 

X n+1 = X„ + v'A^ (A fe (X„) X„ + a k (0)) W„ fc 
fc=i 

+AX n L (X n ) - \ jr X k (X n ) 2 + Ujt\ k (X n ) W% 



k=l \k=l 



+A 6 (0) - £ A fc (X„) a k (0) + E ^ (0) A fe (X„) ll^„ fc 
y k=i i,fc=i 

m(*«)-^E a * W 2 E a * 

fc=l / fe=l 

+A 3 / 2 X„(f;A fe (X„)^ N ) /6 



Vfc=l 



m (*«) - i E A * + M E A * (*«) ^n) E « k (o) ^ fe 

(m \ m 

6(0) -E A " (*») ^ (°) E A " (*«) ^" + ^ (A, X„) , 
fc=l I k=l 

where |ft„ (A,X„)| < |/„| 4 exp (|/„|) | X n + g n | + K (T) A 2 (l + |X„|). This gives 

(x n+1 - x n y - ^ (x„) a + e (*») (wi (n+1) - <„) ^ /& 

< if(T) A 2 (1 + l^j 9 ) (2.16) 

provided that ^ = 1,2,3. 

From (j2~T4|) . pTT5l) and (f27T6j) we obtain (|2~8| . To this end, we can apply the 
classical methodology introduced by Milstein [5^ and Talay [351 [33], or we can directly 
use Theorem 9.1 of [2T] (see also Theorem 14.5.2 of [IS]). □ 
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3. Stable schemes for systems of SDEs. 

3.1. General methodology. Since the coefficients of (jl.ll) are locally Lipschitz 
functions, has a unique continuous strong solution up to an explosion time (see, 
e.g., [29]), which we assume to be +oo a.s. Let b (0) = a 1 (0) = ■ ■ ■ = a m (0) = 0. 
Then, without loss of generality we can suppose Xo ^ a.s. Hence, almost surely, 
X t ^ for all t > 0. 

As we pointed out in Subsection II. 2\ we divide the numerical approximation of 
X t into the computations of \\X t \\ and X t := X t j ||A t ||. Since, almost surely, X t will 
never reach the origin, 

Tj := inf {t > : < — >■,•_►«, oo. (3.1) 



Applying Ito's formula to y ||X t /y TJ || 2 we obtain 



fe=i 



^ ((x a ,b(x.)) + hTZ = i\\° k (x.)\\ 1 f2 {x -' Tk{x 3 s))2 )do 



1**11 2^ \\X„ 



Xt\\ = \\X \\ + f: f {Xs ^ s)) dW s (3-2) 



and so taking limit as j — > oo gives 

* ^ (X s ,6(A s )) + lELilk fc W|| 2 _lf (A s ,a fc (^)) 2 ^ , 
I 11**11 ^ || Xs ||3 ^ 

Moreover, using Ito's formula, together with Tj, we obtain after a long calculation 
that 



X t X f* (b(X s ) I X s b(X s )\ X 



A, || ||A || Jo V II AMI \||AM|' || A, || / || A,.| 



lV T /V Zs * k (X s ) \ 2 / cr k (X s ) a k (X s ) 

A /■*/ X s o k (A s ) \ <r* (A s ) 
"^io \||A S ||' HX.II / ||X S || S 

^ f* (a^Xs) IX^ oMAM \ _A^\ k 

£W V IIAMI \\\x a \\> \\x s \\ J \\x s \\) '• 



ds (3.3) 



I A, 



-ds 



There arc at least two general ways of computing Xx n+1 ■= Xt k+1 / \\Xt, 1+1 ||- For 
example, we can approximate Xt„ +1 by Z n+ i, the numerical solution at time T n+ i of 

Z t = fjnX n + b(Z s )ds + J2 <? k {Z s )dW k , 

JT n fe=1 JT„ 

where fj n and A„ are -measurable random variables such that fj n ~ ||At„|| and 
X n « Xx n '■= Xt„/ \\Xt„ II; as in Sections [T] and [2J for simplicity, we take T n = nA for 
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all n = 0, 1, . . . Then, taking care of round-off errors we set X n+ \ := Z n+ \/ Z n+ \ 

and, in consequence, X n+ \ sa Xt u+1 ■ Alternatively, we can solve numerically (|3.3[) 
which is the method used in Subsections 11.21 and 13.21 In the non-linear case, we can, 
for instance, approximate the right hand side of (|3.3p in each time interval [T„, T n+ i]. 

Now, we provide a manner of handling ||-Xt||. The good performance of the 
numerical method developed in Section [2] motivates us to rewrite (|3.2[) in the form 



\X t \\ = \\X \\+^2 \ 



(X s ,a k (X s )) 
IIA.JI 2 



\XjdW, 



'(X s ,b(X s )) + \\a k (X s )\\ 2 1 ™ (X s ,a k (X s )) 



\X M 



IE 

fc=i 



(3.4) 

\XJ\ds, 



and, further, to approximate on the time interval [T„,T n+ i] by the solution of 

the scalar SDE 



Vt = Vn + Y] \ 



+ 



k=l 
//V 





<7 k ( 


*n)> 






2 









r)s dW s k 



(x n ,b(x n )) + \Y,lli\° k (x f 



\ 



x, 



E 

fe=i 











x n 





(3.5) 



?7 S ds, 



where rj n ~ ||-X^r„|| and X n is an 3r„ -measurable random variable approximating X s 
on [T n ,T n+ i]. Using the explicit solution of (|3.5[) we find 



Xn 



ry„exp 



'(X n ,6(x„)) + iEr=i||^ (*») 



E 

*i=i 



(*») 


2 


m 

E 










fe=i 




a fc ( 


*»)> 






2 





[x n ] 


) 2 






4 



A 



We substitute W# - W$ n by >/AW£, the random variables W$,..., W^, W^,... 
being i.i.d. with symmetric law and variance 1. This leads to 1| fj n +i, where 

fjn+i is given by the recursive formula 



r) n+ i = r\ n exp 



//v 



{Xn,b(x n )) + \Y^ =1 \<j k (x n ) 



x n 



E 

k=l 





a k ( 






x n 





E 

k=l 



(x n , 










2 









AC 



(3.6) 



From (|3.6[) we have that r\ n > for all n S N, in case ?7o > 0. Moreover, if (jl.4p 
and (|1.5[) hold, then we next assert that fj n X n is almost sure exponential stable for 
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all A > provided that X n is a Sr^ -measurable random variable of norm 1. 

Theorem 3.1. Let fj n be as in \3.6\) . Suppose that 6(0) = 0, 770 > 0, and that 
£(770) < 00. Assume that the inequalities and \1.5)) hold. Then 



lim sup — In (t7„) < - A 

71— >oo TLL\ 



(3.7) 



Proof. Deferred to Section f3. 4. II □ 

3.2. Bilinear SDEs. We will construct numerical schemes for bilinear SDEs 
based on p. 21) and p. 31) . This fleshes out the new technique. 

In this subsection, we assume that X t satisfies (|1.2I) . Then, (|3.3p becomes (|1.10[) , 
and so X t is approximated in [T„,T„ + i] by Z t / \[Zt\\, where 



z t = x n + b (y n ) z s ds + E / (° k 



with X n ss Xt and 



B 



B 



Y n ,BY n 



E 

k=l 



k=l ' 



Y n K Y 
1 n 3 u ± ri 



Y n K Y 



Z s dW k 



(3.8) 



Y n ,a K Y n )a K 



a k Y n 



Moreover, (|3.6[) takes the form 

fjn+l = fjn CXp (Y n , BY n ) 



1 m 2 m \ 

sEp^l -E<^>^> 2 A 

k=l k=l / 

m 



(3.9) 



fe=i 



where Y n is an $T n -measurable random variable of norm 1 that approximates X s on 
[T n ,T„ + i]. Thus, we have to specify both Y n and a way of computing Zr n+1 - The 
simplest selection is Y n = X n , together with the use of the Euler-Maruyama method 
to integrate numerically (|3.8|) . This leads to 

SCHEME 3.1. Define recursively X n+ \ = Z n+ i/ L^n+l L where 



Jn+l 



X n +B (X n ) X n A + £ (<j k - (X n , a k X n )) X n VAW k 



(3.10) 



fc=i 



Here Wq, Wq , . . . , W™, W{, . . . are i.i.d. symmetric random variables having variance 
1. Let fj n+ i be given by (|3.9p with Y n = X n . 

Scheme [31] reduces to the method defined by (|1.13p and (|1.15p when W k = . 
Applying Theorem 13.11 we obtain that Scheme 13.11 is almost sure exponential stable 
for all A > under condition (|1.5p . We next establish that it converges weakly with 
order 1. 

Theorem 3.2. Consider T > and f € Cp(M d , M). Let fj n X n be described by 
Scheme \S.1\ with A = T/N , where JV G N. Assume that W k are bounded random 
variables, Xq have finite moments of any order, and that for every g G Cp(R d ,R), 



\Eg(fj X ) - Eg(X )\ < K (1 + E \X \ q ) T/N 



VN € N. 
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Then for all N e 



|E/ (X T ) - Ef (fj N X N ) | < K(T)(l + E\\X \\ q ) T/N. 



(3.11) 



Proof. Deferred to Section |3"X21 □ 

We now focus on (|1.2p with B ill-conditioned. If the matrix B have very differ- 
ent eigenvalues, then we should carefully approximate the term (X S ,BX S ) in (|1.12p . 
This leads us to the problem of finding good candidates for the random variable Y n 
presented in (|3.8[) and (|3.9[) (see Remark 13. 1|) . 

We now develop an alternative strategy to use (|3.8j) and (|3.9[) . which has yielded 
promising results in our numerical experiments. Suppose that X n and p n are $T n - 
measurable random variables such that X n m Xt„, ||-^n|| = 1, and p n w ||-X"t„||- 
Return to 1)1.100 and 1)1.12)1 . and consider the ordinary differential equation 



(3.12) 



Y n (t) = p n X n + / BY n (s) ds Vt g [T n , T n+1 ] 
Since Y n (t) := Y n (t) / \\Y n (t)\\ satisfies 

Y n (t) =Xn + J T (b- (t n (s) , BY„ (s))) Y n (s) ds : 
from (TTTU)) we obtain that for any t *E [T n , X^^i], 

/t tn -t 

* X s ds + Y, (° k - (x s ,a k X s \^ X s dW ; 
k=i 

where [x)j = JXi ^3 (x„, ^X^ /2 - (x„ a k X s ^ a k 



k 

s ' 



a k X s 



/2 . Hence 



X t « Y„ (t) + (Xn) X s ds + Y.J T - (X»> X * dW s 



m „t 

In the spirit of the Euler-cxponcntial schemes given by |24j we consider 



/t m „t 

^>{X n )X s ds + Y. (o- k -(X n ,o- k X n ))X n dW*, (3.13) 



and so X t ~ X t . Approximating the explicit solution of p. 13)1 we deduce that 
X Tn+l ~ U n+ i, where 

U n+1 = cxp (* (X n ) A) (f n / |f„| + (o- k - (X n ,a k X n )) X n JEw£\ . (3.14) 

Here Wq, W§, ■ ■ ■ , W™, W^, . . . are i.i.d. symmetric random variables having variance 
1, and Y n = cxp (BA) X n . This implies Y n / Y n = Y n (T n+1 ). 
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Similarly, using (I1.12[) and Y n (t) 
we can assert that ||X t || « p t for all t £ 



= ||*«|| + f Tn {Y n {s),BY n {s)) Y n (s) 
T n ,T n+ i], where 



ds 



Pt = \\Yn(t)\\ + 



J T ( E {h kl 4 2 / 2 - (Xn,<r k Xn) 2 /2)\ 



p s ds (3.15) 



Y / (X n ,a k X n )p a dW k 
k=i 



Approximating the explicit solution of (|3.15|) we get pt„ +1 ~ Pn+i, with 
Y, 



Pn+1 — Pn 



\k=l 



cxp ( Y, [h kx n\\ /2 - {X n ,a k X n ) 2 ) A + Y( x n,^ k X n )VAWn 

(3.16) 



k=l 



We have derived the following numerical method. 



Scheme 3.2. Let [X n 



n>0 



be given by the recursive formula 



X n +i — U n +i/ \\U n +i\\ 



where U n +\ is as in (|3. 14)) . Moreover, (p n ) n >o is defined recursively by (|3.16p . 

On the other hand, applying the Euler approximation to f|3 . 1 3[) we obtain that 
X Tn+l ~ K+i, where 



V n+ i = %i 



Yn 



* (X n ) X n A + Y(° k - (X n ,a k X n )) X n VAW%. (3.17) 



fe=i 



Then, we can replace U n +i by V n+ i in Scheme [321 This yields 

Scheme 3.3. Define recursively X n+ \ = V n +\ / \\V n +i\\ , where V n +\ is given by 
(|3.17p . Let (p n ) n >o be described by p,16p . 

We next establish the exponential stability of Schemes 13.21 and 13.31 for any step- 
size A, under Condition (|1.7[) applied to (|1.2[) . Here, for brevity reasons, we do not 
address the rate of convergence of Schemes 13.21 and [ 



Theorem 3.3. Consider Scheme 
oo. Suppose that 



or Scheme 



with po > and E (po) < 



A:= sup(x,Bx) + sup V ( ||cr fc a;|| /2 - (x, a k (x)) 2 ) < 0. (3.18) 
IMI=i \\*W=ikTi K ' 



Then 



limsup — ln(p„) < -A 



a.s. 



Proof. Deferred to Section 13. 4.31 □ 

Remark 3.1. Applying the midpoint rule to the solution of (|3.12l) we obtain 



X„ ~ Yn 



cxp (BA/2)X n 
\exp(BA/2)X n \ 



(3.19) 



and so we can approximate (X S ,BX S ) by (Y n ,BY n ) whenever s £ [T„,T n +i]. More- 



over, (X s ,a k X s ) w (Y n ,a k Y n ) and 



a k X s 



\cr k Y n \\ for all s e [T n ,T n+1 ]. 
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Therefore, the choice Y n = Y„ in (|3.8[) and (|3.9[) yields Xj- n+1 ~ ^T„ + il \\Zt„ + i \\ 
and \\Xt„ +1 \\ ~ ^n+i- We can solve (|3.8[) using the Euler- exponential method 
introduced in WM. This gives the recursive scheme X n+ i := V n+ i/ \\V n+ i\\, where 



Vn+i = exp (B (Y n ) A) ix n + ^ (a k - (y„,cr fe F„» A^VAW^J , 

wii/i W™ as in Scheme \3.1\ Finally, we take fj n+ i defined by (|3.9[) w«t/i Y n = Y n . A 
good alternative to V n +i is to apply the backward Euler method to (|3.8j) . 

3.3. Numerical experiments. 

3.3.1. Lyapunov exponents. Wc will illustrate the potential of Scheme 1331 to 
compute the Lyapunov exponents of (jl.2l) (see [4[ [34] for classical theoretical and 
numerical references). To this end, we calcule £ (Xq) := lim t _ i . 00 j In ||X t || , where 

with a, b e R, a > and X 7^ 0. In this well-known test problem [HI [34] , £ (X ) 
docs not depend on the initial condition Xq, and further 



_a + b a -bft* cos (29) exp cos (20)) d9 
C cos (29)) d9 



Following [34] we take a = 1 and b = -2. We choose X = (l/S/2, l/-\/2) . In 
the context of (|1.2j) . there is no loss of generality in assuming ||Xo|| = 1, because 
£ (X ) = I (X Q / \\X q \\) for bilinear SDEs. 

In case B is well-conditioned, we can approximate the Lyapunov exponent £ by 

i N (A) := 1 i n 11^^11 = ^-yJ "v 1 ^"; 1 ^ , 

nk > NA 11 in n\\ NA y \\fj n X n \\ ) 

where fj n X n is given by Scheme 13.11 and ./V is sufficiently large. From (|3.9[) we have 

£ N (A) - ^ En=0 C n (X n ) , With 



^ m m \ m 

fe=i fe=i / fe=i 



£„ (a;) = ^(s,-Bi) + ^ 

For concreteness, we consider W k uniformly distributed on [— -\/3~, v3~] . Then 
£ n+1 (A)=£ n (A)(l 1 ^ < Cn ^ Xn) 



n + 1 J (n + l)A 

In order to evaluate the performance of £ n , we also compute £ by means of the algo- 
rithm £ n introduced in pages 1155 and 1156 of |34) . 

Table 13.11 compares the average of 20 realizations of both £n (A) and £n (A) 
applied to (|3.20[) with a = 10 and cr = 20. We have actually computed In (||X 5 oo|| ) 
since we choose TV = 500/ A. Integrating numerically (|3.21[) wc obtain the reference 
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A 


0.1 


0.05 


0.01 


0.002 


0.001 


0.0001 


a = 10, t= -0.48875 


i N (A) 


30.09701 


26.77603 


-3.71372 


-0.57889 


-0.48503 


-0.48811 


In (A) 


-0.52693 


-0.50851 


-0.49928 


-0.49132 


-0.48930 


-0.48937 


o- = 20, I = -0.49719 


Ma) 


86.56572 


114.26780 


93.37706 


-9.63927 


-2.96739 


-0.48334 


In (A) 


-0.55164 


-0.51945 


-0.50120 


-0.49883 


-0.49798 


-0.49703 



Table 3.1 

Computed values for a final integration time T = ./VA = 500 of the Lyapunov exponent I for 
IMW with a = 1, b = -2, X Q = (l/V2, 1/V2) and different parameters a. 



value i = -0.48875 for a = 10, as well as that I = -0.49719 whenever a = 20. 
Table RP1 shows the very good accuracy of £^. In case a = 10, the relative error 

£n (0-1) — I I \l\ is equal to 0.0781, and so £^ (A) provides a good approximation of 



In (0.05) -£ /\£\ = 0.0448. 



I even if A = 0.1. If we increase the noise to a = 20, then 

In contrast, £jv (0.001) produces the value —2.96739. 

3.3.2. Ill-conditioned matrix B. This subsection addresses the test problem 
(I1.16[) . but with more general drift coefficient B. Indeed, we deal with the numerical 
solution of the SDE 

x - x ° + !l{ b l ° 2 )Wfe <^> 

In order to study cases where B is ill-conditioned, we take 61 = —100 and 62 = 2. We 
also choose a = 4, e = 1 and Xq = (1, 2) T . We have that X t converges exponentially 
fast to since 

-A = sup (x, Bx) + sup I o ^ Ik^ll - y^(^ ; cr k x) 2 
=1 INI =1 V"" k=i 



= max{6i,6 2 } + (e 2 - a 2 ) /2 = -11/2. 

Similar to Subsection 11.21 we compute the mean value of the bounded ran- 
dom variable arctan ^1 + (Xf) 2 j , where X t = (X},X 2 ). In Figure QJ the solid 

line provides the reference values for Earctan^l + (X 2 ) 2 ^, which have been ob- 
tained by sampling 10 s times the Euler-Maruyama scheme E n applied to (|3.22[) with 
step-size A = 2~ 14 w 0.0000610 and Wf n+1)A - replaced by \fK&. Here, 

£o>£o> • ■ • >£cm£i> • • ■ are independent random variables taking values ±1 with prob- 
ability 1/2. Moreover, the circles (resp. stars) represent the estimated values of 

E arctan ^1 + (X 2 ) 2 ^ produced by averaging over 10 6 random realizations of Scheme 
3.21 with W£ = (resp. the backward Eulcr method E n given by (|1.6p with 



^(n+i)A — replaced by s/A£%). Figure 0] shows that the second coordinate 

of E n does not converge to when A = 0.25,0.125, and it goes too fast to for 
A = 0.0625. On the contrary, Scheme l3~2l is very accurate, even with A = 0.25. 

Table 13.21 presents the errors produced at time T = 0.5 by the new numeri- 
cal methods. We assign the weak error e (Y, A) |E/ (Xt) — E/ (Yt/a) | to every 
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£ * # * * * ■* * 




timet 



Fig. 4. Computation of E arctan ^1 + (X^) 2 ^ , where t £ [0,4] and X t solves 113,221 1 with 

bi = — 10, 62 = 2, cr = 4, e = 1 and Xq = (1, 2) T . The true values are plotted with a solid line. The 
circles and stars represent the approximations obtained by Scheme \3.2\ and E n respectively. 

scheme Y with step-size A, where / (xi, x%) = arctan (l + (x 2 ) 2 ). TablcO compares 
estimates of e (Y, A) obtained by sampling 10 6 times the backward Euler method E n> 

Table 3.2 

Absolute errors t(Y,A) = |E/ (Xt) — IE/ (Yt/a) \ involved in the computation of E/ (Xt), 
where X t verifies W1M with 61 = -100, b 2 = 2, o = 4, e = 1 and X = (1, 2) T . Here T = 0.5 and 
/ (3:1,3:2) = arctan ( 1 + (^2) 2 ) ■ 



A 


E 


Scheme [3T2l 


e(Y,A) 
Remark 13.11 


Scheme 13.11 


Scheme 13.31 


1/2 




0.11463 


0.16262 


0.15767 


0.097492 


1/4 


0.62209 


0.033856 


0.039728 


0.15767 


0.033187 


1/8 


0.45494 


0.028848 


0.055127 


0.15767 


0.027073 


1/16 


0.1167 


0.00092922 


0.031429 


0.15767 


0.00001519 


1/32 


0.051448 


0.0034819 


0.02738 


0.15765 


0.0030137 


1/64 


0.022656 


0.0013437 


0.022279 


0.057997 


0.001135 


2- 7 


0.010789 


0.0007307 


0.016977 


0.015613 


0.00065061 


2 -8 


0.0052678 


0.00037132 


0.011165 


0.0061579 


0.00034013 


2 _9 


0.0031195 


0.00034469 


0.0059087 


0.0032785 


0.0003577 


2-10 


0.001551 


0.00017328 


0.0032079 


0.001571 


0.00017914 


2- 11 


0.00045271 


0.00023427 


0.001994 


0.00044898 


0.0002315 


2-12 


0.00031327 


0.000029582 


0.00092728 


0.00030817 


0.000028237 
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Schemes I3.H3.3I and the method introduced in Remark 13.11 We take W k = £™ , and 
the length of all the 99%-confidence intervals are at least of order 10 -3 (see, e.g., [T5]). 
Table I3T21 shows that the values of e ( Scheme 13.21 A) and e (Scheme 13.31 A) arc quite 
similar. We can also see that Schemes 13.21 and 13.31 are very accurate. Moreover, the 
second coordinate given by Scheme 13. II decays too fast to when A > 1/32. Finally, 
the numerical method given by Remark 1 3 . 1 1 has exponentially stable trajectories, but 
tends to slightly more slowly than the 'true' solution. 

3.4. Proofs. 

3.4.1. Proof of Theorem [37TI Using (j3Uj) yields 



In (r) n+ i) =ln(?/ ) + S„ 



E 

3=0 



a k X 



where S n = £™ =1 



(x,-,<r fc (x,-)) 



x, 



E 

k=i 



AW k . Then, (TO) leads to 



{X h a k [X 



X, 



A. 



n + 1 

From dL4l) it follows that 



hi (??„+i ) < ^—ln(?7 ) 



n + 1 



E E 



k=l 



x, 



n + 1 



-S n - AA. 



(3.23) 



X, 



W k 



< KA, 



and so applying a generalized law of large numbers, as in the proof of Theorem 12.1 
we obtain that S n / (n + 1) a.s. By (|3~23)) . 



lim sup ■ 



1 



ln(r?„+i) < -AA 



a.s. 



□ 



n—>oo fl -\- 1 

3.4.2. Proof of Theorem 13.21 Wc first establish that for an arbitrary q>2, 
i 



\v n X n \\ H <K{T)E\fj \ q Vn = 

To shorten notation, we set A fc (X n ) = (X n ,a k X n ) and 



,N. 



(3.24) 



n(X n ) = (X n ,BX n )+J2(\\cr 



k XJ 2 



<X„,a fc X n ) /2. 



fc=i 



Similar to the proof of (|2.14l) . we rewrite as exp (h n ) r] n , and so 

fjn+l = fjn + (cxp {hn) -l-h n )fj n + h n fj n , 
with h n := (X„) - i Efeli A fc (^„) 2 ) A + £™ i (*») VAW* Thus 



?7n+i = ?7o + ^2 (exp (/ife) - 1 - hk) r]k + ^ h k r]k. 

k=0 k=0 
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2 sxp [\n k \) we ODtain 

i 



Using |exp (h k ) — 1 — h k \ < \h k \ exp (\hk\) we obtain 

n n 

\f]n+i\ q <K \fj \ q + K (n+l) q - 1 J2\hk\ 2q \fj k \ q + KA q ' 2 ]T £ X j (X k ) fj k W[ 



k=0 



A (n + If' 1 J2 A 9 



k=0 



i m 

(**)-*£*W 



N 9 . 



k=0 



(3.25) 



For any fceZ +! \\fi (X k ) \\ < K and ||A J ' (A fe ) || < A", because ||X fe || = 1. We also 
have that Wq is a bounded random variable. Then, applying the Burkholdcr-Davis- 
Gundy inequality we deduce from (|3.25[) that 

n 

^\Vn+i\ q < AE| ? 7o| 9 + A(T)A^E|^| 9 

k=0 

for all n = 0, . . . , N — 1. A discrete Gronwall lemma (see, e.g., [3]) now leads to (|3.24[) . 

We proceed to find a truncated asymptotic expansion of fj n+ iX n+ i as A goes to 
0. In what follows, we use the same symbol O (•) for different random functions from 
[0,T] to R or R dxd such that \\0 {x)\\ < A (T)x. Since 

lexp (x)-l + x + x 2 /2 + x 3 /6\ < x 4 exp (|a;|) , 



Vn+1 



Vn (l + L (X n ) - I fj X k (I n ) 2 j A + J2 X k (X n ) W%VZj 

+Vn ^ (Xn) -\Y. XJ (*n) 2 j (ft X " (*») A3/2 



(3.26) 



Vn 

' 2 



6 



J2\ k (X n )W*) A 3 / 2 + ^0(A 2 



\k=l 



\k=l / 

Multiplying the right hand sides of (|3.10[) and (|3.26p yields 

fjn+iZ n+1 = 1 1 + BA + jr o k W k VA J f] n X n + T n A^ 2 f] n X n (3.27) 



fc=i 



+A^T(a fc -A fc (X n ) /2) X k (X n 



k=l 



"7 



l ?y n A„ 



+A ^ {ai - X' (X n ) /2) X k (X n ) W^W k f) n X n + O (A 2 ) fj n X n , 

where T n is a random matrix such that ||r n || < A and E (T n /$T„) = 0; throughout 
the proof, we assume without loss of generality that W^, . . . W™ are 3t„ + i -measurable 
and independent of 3r„ ■ Indeed, 

/ m m \ m 

r « = 5 + ? £ Afe - £ Afe (*») £ AJ (*») ^ 



fc=i 



^ E A " (*0 ^ + U (*») - 1 E Afe (*«) 2 £K- A * (*»)) 

\fc=i / V fc=i / fc=i 
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From (|3.10[) we have 



Z 



n+l 



1 + AJ2 (A" (Xn) Ik^nll ) 1 - (W* 
fc=l ^ 

+A ((a k X n ,a j X n ) - A J (X n ) X k (X n )) W*W£ 



(3.28) 



+ 2A 3/2 ^> ^ ^ ^ _ X k ^ yyk + q ^2 



fc=l 



Hence, there exists Aq > such that 



1- Z, 



n+l 



< A". Using the power series expansion of ir i— > (1 + x) we get 



< 1/2 for all A < Aq, because 

-1/2 



T3 = l + \ 



l-Wz 



1- z 



E-^-fi-IWI 5 



fc-2 



whenever 



l- IU! 



< 1. This, together with (|3.28p . implies that for all A < A , 



i/ p n+1 || = i + f E ( Afe (*«) 2 - h" x 4 2 ) {^) i 



k=i 



(3.29) 



- ^ ((a fe J?„, o*X n ) - X> (X n ) X k (X n )) W k W n 



(3.30) 



-A 3 / 2 ^ (B (X„) X„, <r fc X„ - X k (X n ) X n ) W k + O (A 2 

k=l 

Combining (POT]) with (pT2"9"]) gives 
?7„ + iX n+1 = ^1 + BA + a k W k v^ fj n X n + T n A 3 ^f} n X n 

" l / 2 \ 

+A^(A fc (A > „)a fc - \\a k X n f /i) - 1 ) f) n X n 

k=l ^ ' 

+A J2 0* (X n ) a j - (a k X n , a j X n ) /2) W^W* fj„X„ + (A 2 ) 7j n X„, 

where A < Ao and f n is a random matrix satisfying f„ < K and E (t n /3T„J = 0. 

By flUm , it; is sufficient to prove that (|3TTj) holds for all AT > T/A . Then, from 
now we suppose A < Ao. Looking at (|3.30[) we easily see that ||^„_|_iX n+ i — ^„X„|| < 
K (T) A 1 / 2 ||7?„X n ||, and so 

E(\\fj n+1 X n+1 -fj n X n \\ g /$ Tn ) <K(T)A«/ 2 (l+\\fj n X n \\ q ). (3.31) 
Moreover, (|3.30[) leads to 



( 7?„+iX„ +1 - fj n X n - \ BA + £V ( W T n+ i - w h)) VnX n /3 Tn 
<K(T) A 2 (l + \\fj n X n \\). 



k=l 
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Using again (|3.30[) wc deduce that, up to terms of order O (A 2 ) ||r? n X ra || 9 , the second 
and third moments of fj n +iX n +x — f\ n X n coincide with that of 

m 

Bf] n X n A + <T k fj n X n W*VA, 

k=l 

and then with that of (bA + YJk=i ° k ( W T n+1 - W t„)) VnXn- Therefore, combining 
classical arguments [201 GUI 133] with (|3.24p and (|3.3ip we can assert that (|3.1ip holds 
for all T/N < A (see also Theorem 14.5.2 of [Eg]). □ 

3.4.3. Proof of Theorem [3T3l From (|3~l"6)) we have 



where S n = Y.% a Et!^,^^-)^. Since Y n (t) := cxp (B (t - T n )) X n satis- 
fies (|3.12p . using \\Xj\\ = 1 we get 

ln(||y n (i)||)= f B^-)ds< (t-T n ) sup (x,Bx), 

Jr n \\ y s\\ \\y s \\ \\ x \\=i 

and so In (||exp (BA) X 3 ||) < A sup^n^ (x, Bx). We now apply p,18[) to obtain 

In (pn+i) < In (p ) - (»+ 1) AA + S n , 

As in the proof of Theorems 12.11 and 13.11 a generalized law of large numbers gives 
S n / (n + 1) — > a.s., and the theorem follows. □ 

4. Conclusion. We introduce a new idea to solve numerically multidimensional 
SDEs with multiplicative noise, whose dynamical properties can be meaningfully in- 
fluenced by the noise terms. Applying the new methodology we obtain stable numer- 
ical methods for bilinear systems of SDEs, which have good numerical performances. 
Bilinear SDEs appears in applications such as dynamical studies of the motion of heli- 
copter rotor blades in a turbulent wind. We also develop a promising stable numerical 
scheme for scalar SDEs with multiplicative noise. 
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